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Abstract 



We develop a method to solve the Bogoliubov de Gennes equation for 
superconductors self-consistently, using the recursion method. The method 
allows the pairing interaction to be either local or non-local corresponding to 
s and d-wave superconductivity, respectively. Using this method we examine 
the properties of various S — I — N and D — I — N interfaces. In particular 
we self-consistently calculate the spatially varying density of states and the 
superconducting order parameter. We see that changing the strength of the 
insulating barrier, at the interface, does not, in the case of an s-wave super- 
conductor, dramatically, change the low energy local density of states, in the 
superconducting region near the interface. This is in stark contrast to what 
we see in the case of a D — I — N interface where the local particle density 
of states is changed dramatically. Hence we deduce that in calculating such 
properties as the conductance of S — I — N and D — I — N structures it is far 
more important to carry out a self-consistent calculations in the d-wave case. 
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I. INTRODUCTION. 



The Bogoliubov de Gennes equations and approximations to these equations, for example 
the Eilenberger equations have in the last 10 years been used to great sucsess to examine the 
properties of normal superconducting hybrid structures. Most of this work, in recent times, 
has centered upon solving these equations for the electrical mesoscopic tranport properties of 
such hybrid structures for a recent review of calculations of electrical conductance properties 
of normal superconducting hybrid structures see Lambert and Ramondi [|IJ. Virtually all 
of this work has been centred around the properties of systems which contain conventional 
superconducting regions, ie. the superconducting gap has been isotropic. Also all but a 
few have been non-self-consistent calculations, this however has not inhibited the 

sucsess of such calculations. In more recent times, with the discovery of high-temperature 
superconductivity again the Bogoluibov de Gennes equations have been used to examine 
the properties of hybrid structures, where now the superconducting gap is taken to be non- 
isotropic, usually d-w&ve in symmetry, but mainly so far most of the theoretical analysis of 



simple experiments ||13l-|i5|j have relied upon macroscopic symmetry arguments and not the 
microscopic details of the actual interfaces. However, the microscopic physics at the interface 
can be an important factor in understanding the experimental results. For example, if 
mixing of different order parameter symmetries occurs at the interface (because the interface 
breaks the bulk tetragonal or orthorhombic symmetry) the extent of such mixing can only 
be determined from microscopic calculations. Similarly, suppression of the order parameter 
(either d-wave or s-wave) near an S — N interface can lead to significant local density of 
states within the bulk energy gap, and this can complicate the analysis of single electron 
tunnelling spectra |16| . 

In the past few years there have been a number of microscopic calculations of surfaces and 
interfaces in superconducting systems with a (i-wave order parameter. Most of the theoret- 
ical results have been obtained using tunnelling theory, or Andreev's approximation [|17) -|26| 
in which the tunnelling barrier and the order parameter are not found self-consistently. 
For tunnel junctions these approximations may be adequate, but we show below that self- 
consistency has significant effects for interfaces with direct contact between the constituents. 
Self-consistent properties of interfaces have previously been computed using the Eilenberger 
equations [f2"7]-p2]], which are an approximation to the Bogoliubov de Gennes equation. These 
self-consistent solutions to the Eilenberger equations have shown some interesting effects. 
Also there have now in the last couple of years been a few self-consistent calculations of 
the Bogoliubov de Gennes equations, with non-local interactions, looking at the effect of 
impurity scattering |33|] and superconducting interfaces |34| , [35 . 



In this paper we tackle a simple problem, solving the Bogoliubov de Gennes equations 
self-consistently, in the presence of an s-wave and a <i-wave order parameter with the aim 
of asking, how important is self-consistency? We present a method of how to calculate 
self-consistent properties of superconducting interfaces by directly solving the Bogoliubov 
de Gennes equations [[33] and that certainly in the case of a c?-wave interface it is crucial to 
perform such a calculation to obtain the correct form for the low energy density of states 
in the region of an interface. The method of performing our self-consistent calculations 
will be described in section II. Here we introduce the Bogoliubov de Gennes equation with 
a general interaction, Uij, and demonstrate how this general Hamiltonian can be solved 
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self-consistently using the Recursion Method [|3^,[!7],[3^] . We then present results for bulk 
s-wave and <i-wave systems; looking at both the local particle density of states and in the 
case of d-wave we look at the critical temperature as a function of band filling, we see that 
a d-wave solution, for the interaction we have chosen, is not the only possible solution, for 
our Hamiltonian. 

Then in section III we present results for both S — I — N and D — I — N systems. We 
examine both the form of the superconducting order parameter in the region of the interface, 
for different strengths of insulating barrier, and how the local particle density of states varies 
in the region of the interface. What we see is the sub-gap local particle density of states, 
on the superconducting side of the interface, is not in the case of s-wave superconductivity 
effected by changing the strength of the insultaing barrier, this is in marked contrast to 
what we see in the <i-wave case. Finally in the conclusions we comment on what can be said 
in the light of our calculations. 



II. THE BOGOLIUBOV DE GENNES EQUATION. 

Our starting point is the non-local Hubbard model, which is described by the following 
Hamiltonian 



H— '^^tijC irT Cj IJ [x CjgCiu 4- y y^ j \J^c ia Oj i(T c^_ (J Cj-. a (1) 

ij a ia ija 

where the creation and annihilation operators c\ a and Cj CT , respectivily, create and annihilate 
electrons with spin a in the orbital centred at the lattice point labelled by i, tij is the 
amplitude for hopping from site j to site i, fj, is the chemical potential and Uij is the 
interaction energy of two electrons with opposite spin on sites % and j. To obtain the, 
generalized non-local, Bogoliubov de Gennes Equation from equation ([J) we first make a 
mean field approximation for the pairing field 

Ay U{j (ci a Cj — cr) (2) 

and assuming that the fluctuations about this mean are small we then perform the Bogoli- 
ubov cononical-transformation, enabling us to obtain 

with 

H i6 = ( f ? % ) (4) 

where uf and vf are the particle and hole amplitudes, on site i, associated with an eigenen- 
ergy E n and where Ay is the, possibly non-local, pairing potential or gap function, which 
does not have to be real. 

In the fully self-consistent Bogoliubov de Gennes equation the normal state Hamiltonian 
Hij is given by 
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Hij = {I, j + -UijUij) + (ei - fi)Sij + -Uijrijj (5) 

where /jl is the chemical potential, e« is the normal on site energy of site i and t^, as mentioned 
above is the hopping integral between site i and site j. The interaction terms \Uijfijj and 
"nUijTiij are the Hartree-Fock potentials, which are defined by 

md = E<4q*r) = 2E(«)\ n /(^) - /(^n))- (6) 

Similarly the pairing potentials are defined as 

A i3 = -L U F U (7) 

where the anomalous density is 

Fa = (c iaCj - a ) = £««)*(! - f(E n )) - Ky(u])f(E n ). (8) 

In equations (§) the sum only considers terms _E n up to the condenstate chemical potential 
(/i) as compared to equation (|8|) where the sum is over all terms. 

A solution to the above system of equations will be fully self-consistent provided that 
both the normal (UijUij) and anomalous (A v -) potentials are determined consistently with 
the corresponding densities and F^ via equations (0) and (||). 

Figure 1 illustrates the geometry corresponding to this system of equations, on a tight 
binding lattice. The tight-binding lattice has nearest neighbour hopping interactions (tij), 
as well as a coupling between particle and hole space, via a superconducting order parameter 
(Ay). If the interactions are purely on-site (Uu) attractions then the pairing potential will 
be purely local (A#), corresponding to the dashed line in figure 1. On the other hand when 
the interaction is non-local (Uij, i ^ j) the pairing potential Ay will also be non-local, as 
illustrated, for nearest nieghbour interactions, by the solid lines in figure 1. 



A. The Recursion Method for solving the Bogoliubov de Gennes Equations. 

If we now wish to solve the Bogoliubov de Gennes equations self-consistnetly, we have 
found that the most efficient method of doing this is the Recursion Method |36|. This 
method allows us to calculate the electronic Green's functions 

G aa >(i,j,E) = - H \i a ') ( 9 ) 

where the indices i and j denote sites, while a and a' represent the particle or hole degree 
of freedom on each site. We denote particle degrees of freedom by a = + and hole degrees 
of freedom by a = — . For example G+-(i,j, E) represents the Greens function between the 
particle degree of freedom on site i and the hole degree of freedom on site j. 

To compute the Green's functions (Q) we can closely follow the method described 
by Litak, Miller and Gyorffy P7| , they solved the Bogoliubov de Gennes Equation self- 
consistently, using the recursion method, for a local interaction. Using this method we can 
transform the Hamiltonian to a block tridiagonal form 
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where a n and b n are 2x2 matrices. Given this form for (ia\El — H\ja') and expressing 
the Green's function as 



G aa ,(i,j,E) = (laKEl-Hy'lja'), 
the Greens functions above can be evaluated as a matrix continued fraction so that 



(11) 



G(i, j, E) = \ E1 -a - b\ (El -cn-bl (El - a 2 - b\ (El - a 3 - . . .y 1 63) * 62) ' &i 

(12) 



where 



G(i,j,E) 



G aa (i,i,E) G aa ,(i,j,E) 
(j,i,E) G a , a ,(j,j,E) 



(13) 



Within equations (|T0|) and (|T2| ) we have a formally exact representation of the Green's 
functions. However in general both the tridiagonal representation of the Hamiltonian, and 
the matrix continued fraction (ID) will be infinite. In practice one can only calculate a 
finite number of terms in the continued fraction exactly. In the terminology of the recursion 
method it is necessary to terminate the continued fraction [36- 12 . 

If we were to calculate up to and including a n and b n and then simply set subsequent 
coefficients to zero then the Green's function would have In poles along the real axis. The 
density of states would then correspond to a set of 2n delta functions. Integrated quantities 
such as the the densities and Fy could depend strongly on n, especially since only a few 
of the 2n delta functions would be within the relevant energy range within the BCS cut off, 
E c . In order to obtain accurate results it would be necessary to calculate a large number of 
exact levels, which would be expensive in terms of both computer time and memory. 

As a more efficient alternative we choose to terminate the continued fraction using the 
extrapolation method, as used previously by Litak, Miller and Gyorffy PTfl . We calculate 
the values for a n and b n exactly up to the first m coefficients using the recursion method. 
Then, noting the fact that the elements of the matrices a n and b n vary in a predictable 
manner ||37|| , we extrapolate the elements of the matrices for a further k iterations, where k 
is usually very much greater than m. This enables us to compute the various densities of 
states, and the charge densities n^ and accurately with relatively little computer time 
and memory. 

In terms of the Green's functions G aa '(i, j, E) the pairing and normal Hartree-Fock 
potentials and \Uijnij can be expressed as 



-1 
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and 



Aij = ^ Uij I E + W)- G +-&3, E - ^))(! " V{E))dE (14) 
\Ui?Mi = ^U ij J(G + + (i,j,E + ir ] )-G + + (t,j,E-ir ] ))f(E)dE (15) 



where rj is a small positive number. 

To obtain the above equations we have used the fact that 

■(»?) 




and ^ ( ; » j (16) 

are the eigenvectors of equation @ with eigenvalues E n and — £? n . 

It should be noted that up to this point we have been completely general, so we are not 
just restricted to calculating the Greens functions on the same site, G aa '(i, i, E) or just be- 
tween nearest neighbours, in essence we could use this formalism to calculate a Greens func- 
tion between any two sites. This enables us to also calculate such quantities as G++(i,j, E) 
and G + _(i,j,E), which can then in turn be used to calculate normal reflection and trans- 
mission coefficients, thus enabaling in principle for transport properties or our Hamiltonian 
to be calculated. 



B. Achieving Self-consistency. 



Using the above methods to calculate Ay and Uij-riij we need to achieve a fully self- 
consistent solution. Firstly we make use of any symmetries in the system in order to minimise 
the number of calculations which are necessary. For example on an infinite square lattice 
with no variation in any of the potentials only one independent site needs to be calculated 
since this site can be mapped onto all of the other sites. Secondly, once we have decided 
which sites need to be calculated self-consistently Ay and ny can be calculated for those 
sites, remembering that on a square lattice each site will have four nearest neighbours. This 
implies that in general we will have to calculate nine different Green's functions in order 
to calculate Ay and ray. This can be seen by considering site i in figure 1 and noting 
that we need to calculate the Greens functions shown in table I, depending on whether the 
interaction is purely local, purely non-local, or both local and non-local. Having calculated 
the appropriate Green's functions new values for Ay and n.y can be calculated, which we will 



denote as A; 



(i) 



and n£ . Inserting these into the Hamiltonian and repeating the calculation 



(2) (2) 

of then Green's functions leads to a new set A„V and n)/ and so on. We repeat this iteration 



for all i and j until 



and 



ij 







|A(») 






n (n)| 


n( n ) 



< 0.001 



< 0.001. 



(17) 



(18) 



Since Ay and ny can be complex we need to also check for convergence in their associated 
phases. 
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C. Uniform properties of our Hamiltonian. 



Up to this point we have said nothing about the interaction itself, just that it exists. In 
conventional BCS theory [fE| we can impose a cut-off for our interaction, since we know that 
in this case the interaction is phonon mediated, thus in equations (14) and (16) we can have a 
cut-off for our integrals. However we also wish to consider the properties of high-temperature 
superconductors, where the mechanism for the interaction is a contriversial area, so we state 
at this point that we are not interested in the mechanism for high-temperature superconduc- 



tivity, but we do know, through experiment [44], that the symmetry of the superconducting 



order parameter in high-temperature materials is of a c?-wave nature, therefore, generally 
we wish to consider the propeties of <i-wave superconductors. 

The simplest model to consider is one where we have a nearest neighbour attractive 
interaction, Uij. For such a situation, where the interaction is constant we can calculate the 
critical temperature for a given interaction strength, as a function of band filling, see figure 
2. It is found that there are two possible solutions for the symmetry of the superconducting 
order parameter. In this figure we can see there are two critical temperatures, one for 
extended s-wave solutions, near the edges of the band, and a d-wave solution near the 
centre of the band. It is this (i-wave solution which we will be considering, however it is 
important to note, and will be seen later, that in non-uniform systems although we may be 
at a band filling where only the <i-wave solution exists this does not nessesitate that the 
extended s-wave components to the superconducting order parameter will thence be zero. 

Before proceeding it is instructive to examine the bulk local particle density of states 
(N(E)) for the systems we are concerned with. In figure 3 we have plotted the local particle 
density of states, 

Ni(E) = ^-(G ++ (i, i,E + i V ) - G + + (i, i,E- (19) 

Z7T 

for three different scenarios, no interaction (thin dashed line), a local interaction (thin solid 
line) and a, nearest nieghbour, non-local interaction (thick solid line). What we see from 
this is exactly what one would expect, i.e. for no interaction, we obtain the density of states 
for a regular, nearest nieghbour tight binding lattice, for a local interaction we obtain a BCS 
type gap centred around the Fermi Energy, and finally for a non-local interaction we see a 
typical (i-wave gap, ie. the local particle density of states goes to zero linearly at the Fermi 
energy. Having said that there is no way, at this point to descern that we actually have a 
(i-wave solution, but we define a (i-wave order parameter in the following manner 

I A? } | = \\A in - A ih + A ija - A ij4 1 , (20) 

where ji, j2, js and j± are the nearest nieghbours to i, moving in a clockwise or anticlockwise 
direction, this quantity is then the (i-wave component to the order parameter at i, similarly 
we can also write 

|Ah = \\A in + A ih + A ij3 + AyJ (21) 

for the extended s-wave component to the order parameter. For the results obtained in 
figure 3 lAf ! = 0.38 * t and |A 4 W | = 0. 
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III. RESULTS FOR S-I-N AND D-I-N INTERFACES 



Now we can proceed to perform self-consistent calculations for interfaces. Our model is 
as follows, we consider a system which does not vary in the y coordinate and has a step 
like function in the interaction strength in the x coordinate; ie. for x < the interaction 
is finite; whereas for x > our interaction (be it local or nearest nieghbour non-local) is 
zero. The interface itself is modelled by a line of sites at x — 0, where we change the on- 
site energy of these sites in accordance with the strength of the barrier at the interface we 
wish to model. We then proceed to calculate the Bogoliubov de Gennes equations in the 
self-consistent manner prescribed above. 

First we look at the results obtained for a local interaction (Uu = —1.5*). We see in 
figure 4 the nomalized magnitude of the superconducting order parameter (upper graph) 
and ^U(x)n(x) (lower graph), where in reference to the prevoius equations we have made a 
transformation from sites i to position x. In this figure we see how these two quantities vary 
in the region of the interface (x — 0), for different strengths of the barrier at the interface, 
e(0) = (thin solid line), e(0) = 0.5* (thin dashed line), e(0) = 1.0* (thick solid line) and 
e(0) = 2.0* (thick dashed line). The main point to deduce from this is that the barrier has 
very little effect upon the actual order parameter or the local density. This is in comparison 
to what we see for an equivalent calculation for a non-local interaction, the results of which 
are shown in figure 5, where we have plotted the nomalised <i-wave contribution to the 
superconducting order parameter (upper graph) and the extended s-wave contribution to 
the superconducting order parameter (lower graph), for the same interface configurations as 
used to obtain figure 4. What we see here is that although A^(x) is not seroiusly modified 
in the presence of a barrier, A^ s \x) is. We can see that as we increase the strength of the 
barrier A^(x) in the region of the interface is reduced. 

We now wish to focus our attention on the local particle density of states, near the 
interface. How this quantity changes will be important in deducing whether it is crucial 
to perform self-consistent calculations to obtain quantities which are reliant upon the local 
particle density of states near the interface. As an example we are going to focus on the 
sub-gap density of states and how this changes for different barrier strengths, and whether 
it is then valid to use non-self-consistent calculations in calculating such quantities as the 
sub-gap conductance? 

With this in mind we have in figure 6 plotted the low energy local particle density of 
states for a local interaction (Uu = — 1.5*), at x — —2. We have plotted this quantity 
for different strengths of barrier, e(0) = (thin dashed line), e(0) = 0.5* (thin solid line), 
e(0) = l.Ot (thick dashed line), e(0) = 2.0* (thick solid line) and e(0) = 3.0* (thin dashed- 
dotted line). Again, as in figure 4 there is little change in the quantity we are interested 
in. In fact from this plot we can say that a quantity such as the sub-gap, S — I — N 
conductance depends very little on changes in the density of states and is totally governed 
by the tunnelling conductance of the barrier. 

However again taking the <i-wave case we can see an altogther different story. In figure 
7 we have plotted equivalent local particle density of states as in figure 6, but for a rf-wave 
interface. What we see is that the low energy density of states in the region of the barrier 
is modified significantly, when we change the strength of the barrier, and hence we can 
no longer state that the conductance of an D — I — N interface will be mainly governed 



8 



by the tunnelling conductance of the barrier. This implies that in such calculations, when 
comparing, for example the conductances of different D — I — N interfaces, the details of 
the barrier effect the sub-gap density of states on the superconducting side of the barier and 
hence self-consistency has an important role to play. 

Finally to highlight the point that the low energy local particle density of states in the 
region of an interface is modified significantly, by the strength of the barrier at the interface, 
when the supercunducting region has a d-wave symmetry, we have in figure 8 plotted the 
local particle density of states at x = — 1 (a), x = — 2 (b), x = — 3 (c) and x = —4 (d). 
for two different strengths of barrier e(0) = (thin solid line) and e = 3.0t (thick solid 
line). From this figure we can see the dramatic change in the sub-gap local density of states. 
Particularly in figure 8(a) (nearest the interface) we see that in the region of E = 0, when 
thier is no interface barrier, there is always a significant density of states at the interface, 
as compared to N(0) ~ when we have a finite barrier. 

IV. CONCLUDING REMARKS. 

In this paper we have developed the Recursion Method for calculating the Bogoliubov 
de Gennes equations self-consistently with both local and non-local interactions. Then as a 
simple problem we examined the self-consistent local particle density of states in the region 
of both S — I — N and D — I — N interfaces, in particular we have focused on how the local 
particle density of states is modified as we change the strength of the barrier. What we see 
is that in the case of the S — I — N interface the low energy particle density of states, on 
the superconducting side of the interface, does not change dramatically as we change the 
strength of the insulating barrier. This is in stark contrast to the comparitively dramatic 
change we see when we perform the same calculation for a D — I — N interface. 

These results imply that in the case of <i-wave interfaces the details of the the interface 
are important and hence to get these details correct we should perform a self-consistent 
calculation. If we do a non-self-consistent calculation for the case of the <i-wave interface, 
i.e. just let the order parameter vary in a step like manner at the interface, then we see that 
the low energy density of states on the superconducting side of the barrier is unaltered as 
we change the strength of the insulating barrier, at the interface. This discrepancy, between 
self-consistent and non-self-consistent calculations, arises from the fact that we have a mixing 
of symmetries in the region of the self-consistent d-w&ve interface and when a self-consistent 
calculation is performed we see that this mixing changes as the interface is changed, this 
effect will not be seen in a non-self-consistent calculation. 

These results may well explain why non-self-consistent conductance calculations for hy- 
brid s-wave structures have over the past years been so sucsessful, i.e the problem can simply 
be seen as calculating the tunnel conductance of the barrier, but also raises a warning flag 
that such a bold approach for rf-wave structures may not be as effective, since our calcula- 
tions show that the low energy density of states changes as we change the properties of the 
interface, we believe that the only way to obtain the correct low energy density of states 
behavoir is to perform a self-consistent calculation. 
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FIGURES 



FIG. 1. This is a schematic diagram of a tight binding lattice, with particle and hole degrees 
of freedom, Ajj couples particles on site i to holes on site j. The difference between local and 
non-local pairing is highlighted by the dashed (local pairing) and solid (non-local pairing) lines. 

FIG. 2. The critical temperature vs. band filling for a fixed non-local, nearest nieghbour 
interaction, the dashed line is the extended s-wave solution to the linearized gap equation and the 
solid line is the corresponding d-wave solution. 

FIG. 3. The local particle density of states for a 2-D tight binding lattice, with no interactions 
(Uij = 0) (thin dashed line), with a local interaction (Uu = 1.5*) (thin solid line) and finally with 
a non-local nearest nieghbour interation (Uij = 1.5t) (thick solid line). For this system tij = 1 for 
nearest neighbours only, fi = and = 0. 

FIG. 4. The magnitude of the superconducting order parameter (upper graph), normalised to 
the magnitude at x = -co and l/2U(x)n(x) (lower graph), as a function of x, where we have an 
interface at x = 0, Each of the lines on the upper and lower graphs are for different strengths of 
barrier at x = 0, e(0) = (thin solid line), e(0) = 0.5* (thin dashed line), e(0) = 1.0* (thick solid 
line) and e(0) = 2.0* (thick dashed line). 

FIG. 5. The magnitude of the ci-wave component to the supercunducting order parameter 
nomalized to its value at x = — oo (upper graph) and the s-wave component to the superconducting 
order parameter (lower graph) as a function of x, for different strengths of barrier at x = 0. We 
have e(0) = (thin solid line), e(0) = 0.5* (thin dashed line), e(0) = 1.0* (thick solid line) and 
e(0) = 2.0t (thick dashed line). 

FIG. 6. The local particle density of states, at x = —2, for an attractive local interaction, for 
x < 0. The different lines are again for different strengths of interaction, e(0) = (thin dashed 
line), e(0) = 0.5* (thin solid line), e(0) = 1.0* (thick dashed line), e(0) = 2.0* (thick solid line) and 
e(0) = 3.0t (thin dashed-dotted line). 

FIG. 7. The local particle density of states, at x = —2, for an attractive non-local interaction, 
for x < 0. The different lines are again for different strengths of interaction, e(0) = (thin dashed 
line), e(0) = 0.5* (thin solid line), e(0) = 1.0* (thick dashed line), and e(0) = 2.0* (thick solid line) 
(thin dashed-dotted line). 

FIG. 8. The local particle density of states, at x = — 1 (a), x = —2 (b), x = —3 (c) and x = —4 
(d), for an attractive non-local interaction, for x < 0. The different lines are again for different 
strengths of interaction, e(0) = (thin solid line), and e(0) = 3.0* (thick solid line). 



13 



TABLES 



Interaction Type 


G+-(i,i,E) 


G+±(i,h,E) 


G+±(i,j 2 ,E) 


G+±(i,j 3 ,E) 


G+±(i,u,E) 


Uu 


Y 


N 


N 


N 


N 


Uij(l - Sij) 


N 


Y 


Y 


Y 


Y 


Uu + Uijil-Sij) 


Y 


Y 


Y 


Y 


Y 



TABLE I. This table shows which Greens functions need to be calculated for systems with 
interactions which are local, Uu, non-local, Uij(l — 5 — ij), or both. The site labels correspond to 
the notation of Fig. 1. 
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